﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.IO;
using System.Collections;
using System.Windows.Forms;
namespace vfife_building_algorithm
{
    public static class Solver
    {
        //public Model model;
        private static Hashtable elements;
        
        public static void caculate()
        {
            //- Preprocess
            FileInfo myFile0 = new FileInfo(@"D:\NODE.txt");
            StreamWriter sw0 = myFile0.CreateText();

            FileInfo myFile1 = new FileInfo(@"D:\ELEMENT.txt");
            StreamWriter sw1 = myFile1.CreateText();

            double ET = Model.Caculation_parameters.endTime;    //- Termination time
            double h = Model.Caculation_parameters.h;           //- Size of time step
            double zeta = Model.Caculation_parameters.zeta;     //- damping factor
            double h2 = h * h;

            double c1 = 1 + zeta * h / 2;
            double c2 = 1 - zeta * h / 2;

            elements = Model.Elements;
            Hashtable NodeCollection = new Hashtable();
            double[] EA = new double[elements.Count];
            double[,] ElmtForce = new double[elements.Count,(int)(ET/h)];

            for (int i = 1; i < elements.Count; i++) {
                Element element = (Element)elements[i];
                for(int j=0;j<(int)(ET/h);j++){
                    ElmtForce[i,j]=element.Elementforce.x_exforce[j];
                }
                //Elemnt No i will be marked with an array that corresponds to its ElementForce
                
                Node node2 = element.getNode2();
                Node node1 = element.getNode1();
                if (!NodeCollection.Contains(node1.Id)) {
                    NodeCollection.Add(node1.Id, node1);
                }
                if (!NodeCollection.Contains(node2.Id))
                {
                    NodeCollection.Add(node2.Id, node1);
                }
                EA[i] = element.Material.EA;     //- Area*Young's modulus of member alpha
            }

            double[] mass = new double[NodeCollection.Count];
            double[] x_coord = new double[NodeCollection.Count];
            double[] y_coord = new double[NodeCollection.Count];
            double[] x_velocity = new double[NodeCollection.Count];
            double[] y_velocity = new double[NodeCollection.Count];


            for (int i = 1; i < NodeCollection.Count; i++) {
                mass[i] = ((Node)(NodeCollection[i])).Mass;//- Point mass of point[i]
                x_coord[i] = ((Node)(NodeCollection[i])).Xcoord;
                y_coord[i] = ((Node)(NodeCollection[i])).Ycoord;
                x_velocity[i] = 0;
                y_velocity[i] = 0;
            }
            
            //Node(i)'s initial coodinate.
            

            
            
            

            
            
            

            
            

            
            
            

            
            
            
            ////- Initial internal force
            //double f1x = 0.0;  //- x-internal force of pt 1
            //double f1y = 0.0;  //- y-internal force of pt 1
            //double f2x = 0.0;  //- x-internal force of pt 2
            //double f2y = 0.0;  //- y-internal force of pt 2
            //double f3x = 0.0;  //- x-internal force of pt 3
            //double f3y = 0.0;  //- y-internal force of pt 3

            ////- Initial external force
            //double curve_value = ramp_step(rtime, time);    //- ramp-step curve w.r.t. time
            //double p1x = 0.0;                               //- x-external force of pt 1
            //double p1y = -Pa * curve_value;                    //- y-external force of pt 1
            //double p2x = Pb * curve_value;                    //- x-external force of pt 2
            //double p2y = -Pa * curve_value;                    //- y-external force of pt 2
            //double p3x = Pb * curve_value;                    //- x-external force of pt 3
            //double p3y = 0.0;                               //- y-external force of pt 3

            ////- Evaluate x_1
            //double x2n_1 = x2n - (h + 0.5 * zeta * h2) * v2xn + 0.5 * h2 * (p2x + f2x) / mass2;
            //double y2n_1 = y2n - (h + 0.5 * zeta * h2) * v2yn + 0.5 * h2 * (p2y + f2y) / mass2;

            //double x2n1;
            //double y2n1;
            //double d2xn1;
            //double d2yn1;
            //double lt;
            //double l0;
            //double fa;
            //double exa;
            //double eya;
            //double fb;
            //double exb;
            //double eyb;

            ////- Timexae integration

            //for (time = 0.0f; time < ET; time = time + h)
            //{
            //    //- Central difference stencile
            //    x2n1 = (2.0 * x2n - c2 * x2n_1 + h2 * (p2x + f2x) / mass2) / c1;
            //    y2n1 = (2.0 * y2n - c2 * y2n_1 + h2 * (p2y + f2y) / mass2) / c1;

            //    double n1 = 1;
            //    double n2 = 2;
            //    double n3 = 3;
            //    double e1 = 1;
            //    double e2 = 2;
            //    double zero = 0;


            //    sw0.WriteLine("  " + time.ToString("e16") + "  " + n1.ToString("e16") + "  " + x1.ToString("e16") + "  " + y1.ToString("e16") + "  " + zero.ToString("e16"));
            //    if (x2n1 >= 0 && y2n1 >= 0)
            //    {
            //        sw0.WriteLine("  " + time.ToString("e16") + "  " + n2.ToString("e16") + "  " + x2n1.ToString("e16") + "  " + y2n1.ToString("e16") + "  " + zero.ToString("e16"));
            //    }
            //    else if (x2n1 < 0 && y2n1 >= 0)
            //    {
            //        sw0.WriteLine("  " + time.ToString("e16") + "  " + n2.ToString("e16") + " " + x2n1.ToString("e16") + "  " + y2n1.ToString("e16") + "  " + zero.ToString("e16"));
            //    }
            //    else if (x2n1 >= 0 && y2n1 < 0)
            //    {
            //        sw0.WriteLine("  " + time.ToString("e16") + "  " + n2.ToString("e16") + "  " + x2n1.ToString("e16") + " " + y2n1.ToString("e16") + "  " + zero.ToString("e16"));
            //    }
            //    else
            //    {
            //        sw0.WriteLine("  " + time.ToString("e16") + "  " + n2.ToString("e16") + " " + x2n1.ToString("e16") + " " + y2n1.ToString("e16") + "  " + zero.ToString("e16"));
            //    }
            //    sw0.WriteLine("  " + time.ToString("e16") + "  " + n3.ToString("e16") + "  " + x3.ToString("e16") + "  " + y3.ToString("e16") + "  " + zero.ToString("e16"));

            //    sw1.WriteLine("  " + time.ToString("e16") + "  " + e1.ToString("e16") + "  " + n1.ToString("e16") + "  " + n2.ToString("e16"));
            //    sw1.WriteLine("  " + time.ToString("e16") + "  " + e2.ToString("e16") + "  " + n2.ToString("e16") + "  " + n3.ToString("e16"));



            //    d2xn1 = x2n1 - x2;
            //    d2yn1 = y2n1 - y2;

            //    //- update time integration data
            //    hc = hc + h;
            //    x2n_1 = x2n;
            //    y2n_1 = y2n;
            //    x2n = x2n1;
            //    y2n = y2n1;

            //    //- Internal forces
            //    //- Elem. [1] (1,2)
            //    lt = Math.Sqrt(Math.Pow((x2n - x1), 2.0) + Math.Pow((y2n - y1), 2));  //- element length @ time = t
            //    l0 = Math.Sqrt(Math.Pow((x2 - x1), 2) + Math.Pow((y2 - y1), 2));    //- element length @ time = 0
            //    fa = EAa * (lt - l0) / l0;                     //- element force
            //    exa = (x2 - x1) / l0;                      //- element orientation @ time 0
            //    eya = (y2 - y1) / l0;                      //-

            //    //- Elem. [2] (2,3) (using the code of Elem. [1] as a templet)
            //    lt = Math.Sqrt(Math.Pow((x3 - x2n), 2) + Math.Pow((y3 - y2n), 2));
            //    l0 = Math.Sqrt(Math.Pow((x3 - x2), 2) + Math.Pow((y3 - y2), 2));
            //    fb = EAb * (lt - l0) / l0;
            //    exb = (x3 - x2) / l0;
            //    eyb = (y3 - y2) / l0;

            //    //- Assemble element internal force to pt
            //    f1x = -(-fa * exa);
            //    f1y = -(-fa * eya);
            //    f2x = -(fa * exa - fb * exb);
            //    f2y = -(fa * eya - fb * eyb);
            //    f3x = -(fb * exb);
            //    f3y = -(fb * eyb);

            //    //- External force @ time + h
            //    curve_value = ramp_step(rtime, time + h); //- ramp-step curve w.r.t. time
            //    p1x = 0.0;                                //- x-external force of pt 1
            //    p1y = -Pa * curve_value;                     //- y-external force of pt 1
            //    p2x = Pb * curve_value;                     //- x-external force of pt 2
            //    p2y = -Pa * curve_value;                     //- y-external force of pt 2
            //    p3x = Pb * curve_value;                     //- x-external force of pt 3
            //    p3y = 0.0;                                //- y-external force of pt 3
            //}
            //sw0.Close();
            //sw1.Close();

        }
    }
}
